
%% initialize
clear all;
close all;

rng(328120)

Nsim = 500;

loadstarting = 0;
loadestimates = 0;

%% read in data

% choices

rawdata = csvread('../../01_data/02_basedata/estimationdata/All_T_J/All_T_J_choices_buff_0_0_lf_principalExtra_ach.csv',1,0);

IJjtidx         = rawdata(:,1);
IJsidx          = rawdata(:,2);
IJpositive      = rawdata(:,3);
IJVAmadev       = rawdata(:,4);
IJVAmean        = rawdata(:,6);

IJX             = [IJVAmadev+IJVAmean,rawdata(:,5)+rawdata(:,7),rawdata(:,8:end)];

clearvars rawdata;

rawdata = csvread('../../01_data/02_basedata/estimationdata/All_T_J/All_T_J_choices_buff_0_0_lf_principalExtra_cf_ach.csv',1,0);

IJitidx         = rawdata(:,1);
IJjtidx_teach   = rawdata(:,2);
IJjtidx_princ   = rawdata(:,3);
IJapp_year_cf   = rawdata(:,4);
IJVAmadev_cf    = rawdata(:,5);
IJVAmean_cf     = rawdata(:,7);

IJX_princ       = [IJVAmadev_cf+IJVAmean_cf,rawdata(:,6)+rawdata(:,8),rawdata(:,9:end)];


clearvars rawdata;

rawdata = csvread('../../01_data/02_basedata/estimationdata/All_T_J/All_T_J_choices_buff_0_0_lf_principal_cf_n_ach.csv',1,0);

IJN_m1_DISAD     = rawdata(:,1);
IJN_m2_DISAD     = rawdata(:,2);
IJN_m1_race     = rawdata(:,3);
IJN_m2_race     = rawdata(:,4);
IJN_m1_ach      = rawdata(:,5);
IJN_m2_ach      = rawdata(:,6);
IJVA_m1_DISAD     = rawdata(:,7);
IJVA_m2_DISAD     = rawdata(:,8);
IJVA_m1_race     = rawdata(:,9);
IJVA_m2_race     = rawdata(:,10);
IJVA_m1_ach      = rawdata(:,11);
IJVA_m2_ach      = rawdata(:,12);
exper_potential = rawdata(:,13);

clearvars rawdata;



%% prepare data for estimation

ds = max(IJsidx);
d2 = max(IJjtidx);



%% draw random effects


Xdraws = zeros(size(IJjtidx,1),1,Nsim);

drawss = normrnd(0,1,[d2,Nsim]);

Xdraws = zeros(size(IJjtidx,1),1,Nsim);
Xdraws(:,1,:) = drawss(IJjtidx,:);

%% estimate mixed logit


y_logistic = IJpositive;


Xmat = [ones(length(IJpositive),1),IJX];

f =  @(x)func_mixed_logit(x,y_logistic,Xmat,Xdraws,IJjtidx,0);


    betastart = zeros(size(Xmat,2),1);
    sigmastart  = ones(size(Xdraws,2),1);
    
    betastart(1) = -1;
 
    thetastart = [betastart;sigmastart];




optionsgrad = optimset('MaxFunEvals',10000000,'MaxIter',10000000,'TolX',1E-5,...
    'Display','Iter','TolFun',1E-7,'Algorithm','trust-region','GradObj','on');


    tic
thetahat = fminunc(f,thetastart,optionsgrad);
toc
f =  @(x)func_mixed_logit(x,y_logistic,Xmat,Xdraws,IJjtidx,1);

[obj,grad,se] =  f(thetahat);

 save '../../01_data/05_temp/mixed_logit_principal_totVA_ach'
 

